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Abstract 

An approach for solving the compressible Euler and Navier-Stokes equations upon meshes composed of nearly arbitrary 
polyhedra is described Each polyhedron is constructed from an arbitrary number of triangular and quadrilateral face 
elements , allowing the unified treatment of tetrahedral, prismatic, pyramidal and hexahedral cells, as well the general cut 
cells produced by Cartesian mesh approaches . The basics behind the numerical approach and the resulting data struc- 
tures are described. The accuracy of the mixed volume grid approach is assessed by performing a grid refinement study 
upon a series of hexahedral, tetrahedral, prismatic and Cartesian meshes for an analytic inviscid problem . A series of 
laminar validation cases are made, comparing the results upon differing grid topologies to each other, to theory and 
experimental data . A computation upon a prismatic/tetrahedral mesh is made simulating the laminar flow over a wall/cyl- 
inder combination. 


I Introduction 

Unstructured grids are rapidly becoming more useful for 
the simulation of inviscid flows in complex geometries. 
The promise of easing the burden of grid generation for 
complex geometries is being met. By exploiting certain 
geometric properties of tetrahedra and convex unit aspect 
ratio hexahedra (Cartesian cells), efficient methods can be 
found that fill the volume of the domain, with some user 
intervention still needed to provide guidance upon ceU 
size and possibly stretching directions. Although the vol- 
ume grid generation can be relatively automated, the sur- 
face discretization of complex geometries is still a non- 
trivial task. There are presently two separate camps of 
unstructured volume grid generation: tetrahedral and Car- 
tesian based. Tetrahedral based mesh generation 
approaches currently being investigated can be grouped 
into advancing-front [1], advancing-layer [2], and point 
insertion [3] methods. Cartesian mesh generation is a rel- 
atively newer approach, which uses a recursive subdivi- 
sion of convex, unit-aspect ratio Cartesian cells, and 
creates (possibly) non-convex polyhedra near boundaries 
[4,5,6]. 

The use of tetrahedral elements can provide efficient cell- 
centered and vertex-based schemes. For a cell-centered 
approach, where the conservation volumes are the tetra- 
hedra themselves, the fixed number of faces and vertices 
of the control volume results in a simpler flow solver. For 
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vertex-based schemes, where the control volume is the 
dual mesh, elegant formulations can result using a conser- 
vative, finite-element framework. By using linear finite- 
elements and exploiting certain geometric properties of 
the tetrahedra, efficient edge-based schemes can be formu- 
lated. The use of tetrahedral grids does, though, have its 
drawbacks. One stems from requiring that the surface dis- 
cretization match faces in the volume grid exactly, which 
makes the surface discretization a controlling part of the 
quality of the volume grid generation. In addition, the vol- 
ume grids generated are irregular in the sense that the ori- 
entation of the faces of the volumes do not typically 
follow a preferred direction. 

Cartesian based approaches attempt to overcome these 
two problems by filling the volume with regularly ori- 
ented, nearly isotropic cells, that become general polyhe- 
dra near the boundaries, where these boundary cells have 
been cut from the Cartesian/boundary intersections. This 
has essentially sacrificed grid smoothness at the boundary 
for grid smoothness over the larger portion of the volume. 
Other benefits of the Cartesian approach can be traced to 
taking advantage of the geometric regularity of the un-cut 
cells, and other implementation specific benefits resulting 
from the hierarchy of the grid from the grid generation 
process. The lack of grid smoothness along the boundaries 
can cause problems for both inviscid and viscous calcula- 
tions, and the resulting solvers are slightly more compli- 
cated than those based upon tetrahedral cells. These 
drawbacks aside, the Cartesian approach is proving to be a 
very successful method for computing inviscid flows 
about complex geometries. 



Both tetrahedral and Cartesian strategies are lacking 
when computing viscous flows. The current viscous flux 
formulations dictate that smoothly stretched, nearly 
orthogonal grids are needed to provide robust and accu- 
rate predictions of viscous flows. The requirement of 
grid smoothness rises to extreme importance, since non- 
smoothness has a direct effect upon needed derivative 
quantities at walls, such as skin friction and heat trans- 
fer, which are typically the quantities desired from such 
an analysis. Grid smoothness also has a direct effect 
upon convergence behavior, since the typical flux func- 
tions in use today will produce non-positive stencils if 
certain geometric qualities of the mesh are not met [3,4]. 
To predict skin friction and heat transfer properly in tur- 
bulent flows, high resolution is needed normal to the 
wall dictating large numbers of cells. In addition, from 
an efficiency standpoint, grid stretching is typically 
needed in only a single direction, normal to the stream 
surface, and is not needed along it. By construction, 
Cartesian based methods do not allow for anisotropy of 
the mesh, while the efficiency of using highly stretched 
tetrahedral cells is suspect. 

A means that is proposed to alleviate these deficiencies 
is currently being called a prismatic grid approach. In 
this case, bounding surfaces are triangulated, and this 
bounding triangulation is extruded away from the sur- 
face, creating layers of cells that are smoothly stretched 
in a surface normal direction. Within the layers, the 
desired smoothness and near-orthogonality is retained. 
These prismatic cells are typically grown out a distance 
from the surface, then a volume mesh generation strat- 
egy is used to fill the void. Examples of this approach 
are shown by Melton et al. [7] for the Euler equations, 
where a Cartesian grid was used to fill the void, and a 
hyperbolic-like approach was used to generate the pris- 
matic layers. Karman [6] used a similar Cartesian/pris- 
matic approach for the Euler and Navier-Stokes 
equations, where a more algebraic approach was used 
for the prisms. Connell et al. [8,9] used a surface discret- 
ization coupled with a CAD-based surface description, 
from which an algebraic approach extruded the prisms, 
and an advancing front mesh generator filled the void 
with tetrahedra. Kallinderis et al. [10] used a similar 
approach, but did not create the surface description from 
a CAD basis. By exploiting the semi-structured nature 
of the prismatic portions of the grid, Parthasarathy et al. 
[11] have proposed an efficient strategy to solve the 
Euler and Navier-Stokes equations. Some obvious 
drawbacks of the prismatic approach, in general, still 
require some work to resolve. For instance, the bound- 
ary surface discretization will control the smoothness of 
the grid near the wall, and care must be taken to ensure 
smoothness at the prism/volume grid interface. Regard- 


less, current examples of this approach show tremen- 
dous potential, where it is hoped to alleviate many of the 
problems unstructured grid approaches encounter for 
computing high Reynolds number, turbulent flows. 

A common thread to computing flows upon these 
classes of grids is that the flow solver must handle both 
tetrahedral, pentahedral (prismatic and pyramid) and 
hexahedral cells. Additional capability to handle adap- 
tive mesh refinement, “hanging nodes”, Cartesian gen- 
erated grids with their cut cells, the extrusion of 
quadrilateral cells into hexahedra, or, perhaps, extrusion 
of other surface polygons would also be desired. In gen- 
eral, this type of solver must be able to solve the conser- 
vation laws upon general, non-simplicial conservation 
volumes. 

The use of edge-based data structures have been pro- 
posed to solve the Euler/Navier-Stokes equations on 
mixed-element meshes by Mavriplis et al. [12]. In this 
case, a convincing argument is made for the use of 
mixed-element meshes, and computations using differ- 
ing element types for the same meshes are performed, 
rather impressively. In [12] the edge-based formulation 
is also used for the discretization of the viscous terms, 
which analysis shows to be inaccurate on non-simplicial 
meshes. An argument is made that relates this discreti- 
zation to a thin-layer like formulation, so for certain 
flows, the results might be adequate, but in general, a 
different formulation for the viscous terms might be 
desirable. This will undoubtedly not be solely edge- 
based, but a careful implementation should not detract 
too much from the approach. 

The approach presented here solves the Euler and 
Navier-Stokes equations using a cell-centered, finite- 
volume scheme upon control volumes of nearly arbi- 
trary polyhedra constructed from triangular and quadri- 
lateral faces. The four basic cell types of tetrahedra, 
prisms, pyramids and hexahedra are a subset of this, 
plus the approach can compute flows upon Cartesian 
generated grids, and grids where cell refinement has 
introduced “hanging nodes”. It is certain that by restrict- 
ing the mesh to be comprised of simpler polyhedra, a 
simpler flow solver results. The approach here is based 
upon the premise that by placing less restriction upon 
the topology of the mesh, an overall faster turnaround 
time will result. The additional computational complexi- 
ties associated with this approach are tractable with well 
thought out data structures and algorithms. One noted 
difference from this approach and standard cell-centered 
methods is that both cell-averaged data and data at the 
vertices of the control volumes are used. 

The outline of this paper is as follows. The basic data 
structures used for the approach are explained, then the 
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approaches used to solve the conservation laws are 
shown. The issues regarding vectorization are 
addressed, namely the coloring of the different opera- 
tions upon the basic data types. The accuracy of the 
approach is assessed using an analytic solution to the 
Euler equations for the four basic cell types mentioned 
above, and Cartesian generated grids. The laminar flow 
over a selection of problems are then made, and compar- 
ison is made to theory (flat plate) and experiment 
(developing duct flow). To demonstrate the capability of 
computing upon prismatic/tetrahedral meshes the 
approach is used to compute the laminar flow about a 
right circular cylinder upon a flat plate using a grid gen- 
erated by Connell et al. [8]. 

IT Data and Data Structures 

The compressible Euler and Navier-Stokes equations 
are solved in three-dimensions in a cell-centered, finite- 
volume format upon a mesh of polyhedra where each 
polyhedron is created from an arbitrary number of trian- 
gular and quadrilateral face elements. This particular 
finite- volume approach uses data at both cell centers and 
mesh vertices, which implicitly uses all nearest neighbor 
cell data, without having to store its connectivity. Since 
the governing equations are being solved in conserva- 
tion law form using a cell-centered scheme, the data 
structures used in the code are designed for such an 
approach. Three separate data entities are identified that 
make up the mesh and are needed to perform the calcu- 
lations: vertices, faces, and cells. These form a hierar- 
chical-like relationship with each other, since faces are 
constructed from vertices and cells are constructed from 
faces. Edges can be obtained from faces/cells, but since 
they are not used directly in this cell-centered scheme, a 
separate data entity for them is not maintained. The sep- 
arate operations needed in the solution of the conserva- 
tion laws are cast as operations upon these data types. A 
collection of vertices/faces/cells which make up a por- 
tion (either complete or by some geometric decomposi- 
tion) of the entire domain are further grouped together 
into data entities called grids. Due to the hierarchy of the 
data entities, much information can be obtained directly, 
such as cell-face-cell connectivity. To clarify these data 
entities, the following briefly describes each entity, the 
data it contains, and how it is stored and maintained. 

ILa Grids 

Grids are composed of lists of cells/faces/vertices upon 
which the computations are performed, as well as lists 
of boundary faces and boundary vertices. All loop color- 
ing information is also stored here. It is intended to 
eventually perform multigrid acceleration, where each 
grid in the sequence will be constructed via an agglom- 


eration approach. This approach lends well to a hierar- 
chy of grids, corresponding to the agglomeration of the 
parent grid and the construction of further grids in the 
sequence via solution adaptive mesh refinement. This 
will be represented hierarchically in the grid data struc- 
ture also, so that each grid will have pointers to its par- 
ent and children. Although multigrid acceleration is not 
employed at this stage of the solver development, the 
use of the grid Hata entity should aid in its implementa- 
tion. The grid data type might also be useful for parallel 
computations, by holding the spatially decomposed 
data. 

n.b Cells 

Cells define the conservation volumes upon which the 
conservation laws are solved and provide a place to 
store cell-averaged data, flux balances and other cell- 
based quantities. Each cell is comprised of an arbitrary 
number of faces, where a list of pointers for each face is 
maintained in the cell. Cells are accessed by a list of 
pointers to the cell data structures. For vectorization of 
cell-to-vertex scatter operations, cells are grouped into 
like vertex number groups and colored (see Section IV). 

ILc Faces 

Fluxes are integrated across faces, and the result of the 
integrations are scattered to the cells which are logically 
left and right of the face where the logical orientation is 
determined by the face normal vector. Faces contain 
pointers to the cells that are logically left and right of the 
faces. Each face contains a list of pointers to globally 
unique vertices defining the geometry of the face. For 
practicality, the faces must be either triangular or quad- 
rilateral. The needed geometric data for the flux integra- 
tion (Gauss points, area vector) are stored in the face 
data structure, but may also be computed from the face 
vertex data. Cell faces are constructed from the posi- 
tively (outward pointing normal) ordered vertex lists of 
the input cell definition. To create a global list of unique 
faces, an octree-based, fixed bucket size searching pro- 
cedure is used to match already created faces. Faces are 
sorted according to Gauss point location in space, and 
are matched by their ordered vertex lists. The tree auto- 
matically sizes itself during the creation phase, and 
prunes itself to zero length when it is not needed. For a 
multi-block grid or a grid with different grid types in 
each domain, there is no need to supply inter-block or 
inter-grid connectivity data, since this face matching 
procedure will automatically ensure the proper face 
matching across inter-block or cell-refinement bound- 
aries. Interior and boundary faces are maintained and 
processed in separate lists. The face data is not only 
used for flux evaluation, but are also used in the upwind, 
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inviscid reconstruction procedure explained in Section 
HI.b. For vectorization purposes, the face loops axe also 
colored (Section IV). For boundary condition applica- 
tion, ghost cells are created for all boundary faces, and 
the data used in the cell-to-vertex scatters. Boundary 
condition faces of the same boundary condition type are 
grouped together, and these are also colored. 

Il.d Vertices 

The vertex data structures hold the spatial coordinates of 
the vertices of the mesh, and provides a list of pointers 
to cells which have edges of faces that are subtended by 
the vertex. This provides a means of obtaining the solu- 
tion at the vertices from the cell-centered data, which is 
needed to compute the upwind, inviscid reconstruction, 
to form the viscous fluxes, and to plot the solution. Ver- 
tex data is obtained in a nominally linearity-preserving 
manne r, as shown in Section IILa. As in the face data, 
upon in put, a self-expanding, bucket searching octree 
procedure is used to provide unique vertex data, where 
sorting and matching is made according to the spatial 
location of the vertex. This makes multi-block and 
multi-grid type data definition easier, since the burden 
of vertex matching is taken by the octree, and not the 
grid generation. 

ILe Tnniit Data Types 

The flexibility of the conservation volume construction 
is evident by the various standard grid data types that 
can be read in by the code. The solver is presently con- 
figured to read in 5 types of grid data; PLOT3D data 
[13], VGRID data (see [2] and others), an input format 
corresponding to the prismatic/tetrahedral grid genera- 
tion system described in [8,9,14], and a format corre- 
sponding to the tetrahedral generator of the FEUSA 
system[15]. Another gnd type is also available, termed 
here as the MVG (Mixed Volume Grid) type, which 
defines each cell as being constructed of an arbitrary 
number of triangular and quadrilateral faces. Cartesian 
generated grids are input as the MVG data type, since 
the cut cells generated by the Cartesian grid generator 
produces polyhedra that are not of the four types listed 
above. All of the 5 specific data types can be translated 
into the MVG format 

TTT Solution of the Conservation Laws 

The Euler and Navier-Stokes equations cast in conserva- 
tion law form are solved in a cell-centered, finite-vol- 
ume format upon the polyhedral control volumes. Both 
upwind and central-difference approximations of the 
convective fluxes are used, and a directed gradients pro- 
cedure is used for the viscous fluxes. A simple three- 


stage explicit scheme with a spatially varying time step, 
based upon both hyperbolic and parabolic time step con- 
straints, is used to update the solution. Contrary to most 
cell-centered approaches, the solution procedure here 
relies upon both vertex and cell-averaged data. The use 
of both vertex and cell-averaged data is also used in the 
USM3D series of unstructured mesh solvers developed 
by F rink [16], and has inspired some of the generaliza- 
tions to mixed-volume grids, shown here. 

m a Vertex Data Interpolation 

The data at the vertices of the mesh is obtained from the 
cell-averaged data by a linearity-preserving interpola- 
tion procedure. This interpolation procedure is similar to 
that presented by Holmes and Connell in [17], where it 
is termed a linearity-preserving Laplacian weighting. By 
considering an interpolation formula that interpolates 
the solution at the j-th vertex from the n cells surround- 
ing it as 

f, - 2V» (1 > 

n 

where the 05 are found from the weights of a pseudo- 
n 

Laplacian operator, 

Lif) = I > n (f n ~fj) ( 2 > 

n 

as 


n 

Linearity preservation is guaranteed by requiring (2) 
satisfy 

L(x) = L(y) = L(z) = 0. (4) 

By expanding the weights about unity in terms of linear 
basis functions, as 

% - i+VV^ + V?,rV + V*n-y< 4 5 > 

a 3x3 linear system is found for the X. which is inverted 
beforehand. This process is purely geometric, and there- 
fore is precomputed for a given mesh, and only requires 
a simple cell-to-vertex scattering procedure. By provid- 
ing higher-order constraints in (4) and expanding with 
higher-order basis functions in (5), quadratic-preserving 
reconstructions can also be found [4] but are not used 
here, since only linearity-preserving cell- and face- 
based reconstructions are used. In practice, the weights 
(5) are restricted 0 < a> n < 2 . 
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Ill.b Unwind Flux Construction 

The upwind scheme follows standard practice used for 
unstructured grids: A reconstruction procedure is used 
to reconstruct the solution locally within each cell, 
which is then followed by an upwind flux construction 
at the cell faces. The reconstruction within each cell is 
used to provide states to an approximate Riemann 
solver, which is used to compute the fluxes. The fluxes 
are then scattered to the logical left and right cells of the 
face. 

nib.l Reconstruction 

The reconstruction of the solution within each conserva- 
tion volume is found by performing a surface integral 
over the conservation volume itself. Since the data at the 
vertices of the conservation volume are found in a lin- 
earity preserving manner, a second-order, Gaussian 
quadrature guarantees that the reconstructed solution is 
also linearity preserving. For a general control volume, 
consider finding the gradient in the reconstruction 

“ = a+ tr. ( *i _ *« } (6) 

by applying a divergence integral over the volume 

= <j)HA^d(3V) * ^ ^ 

V 1 BV 

By replacing the surface integral in (7) with a numerical 
quadrature, the gradient can be found by a face based 
scatter operation. The surface integral is replaced by a 
single point quadrature so that the gradient is found to 
be 


— = - y « 

dx v 2 , u 


faces 


G N i 


( 8 ) 


where N. is the i-th Cartesian component of the out- 
ward facing (non-unit) normal, and Uq is found by 
evaluating a linear or bi-linear expansion in the face at 
the face Gauss point for three- and four-vertexed faces, 
respectively. The cell volumes and centroids are found 
beforehand by using a third-order Gaussian quadrature 
of another application of the divergence integral. For the 
general flux G . , application of the divergence theorem 
replaces the volume quadrature with a surface integral. 
That is 


fcKr. r 

\d7 dV= f G «V( av ) ( 9 > 

V * dV 


To obtain the cell volume and centroid locations the flux 


on the right hand side of (9) is taken to be 


(x, 0, 0) 

for V 

G = 

&*•) 

for xV 


(xy, 0, 0) 

for>V 


_(xz, 0, 0) 

for zV 


The reconstruction (8) preserves linear data, and for a 
general mesh implicitly includes all local cell data in the 
construction of the gradient through the vertex data 
interpolation. Since the faces are colored for vectoriza- 
tion of the integration of the inviscid and viscous fluxes, 
this procedure follows closely the overall face-based 
solution method. 

It is instructive to note that when applied to a mesh of 
tetrahedra, this procedure results in the same formula 
found by Frink [16], also pointed out by Mitchell [18]. 


IIIb.2 Upwind Flux Construction 

Three popular upwind flux functions are available: 
Roe’s FDS [19], Van Leer’s FVS [20] and Liou’s 
AUSM+ [21]. The reconstruction provides the left and 
right input states at the faces, which then use the approx- 
imate Riemann solvers noted above to compute the flux. 
The upwind fluxes are then scattered to the cells. 

IILc Central Differencing with Explicit. 
Scalar Dissipation 

A conservative formulation corresponding to centrally 
differenced fluxes with blended second- and fourth- 
order dissipation is also available in the solver. Follow- 
ing [22], at a given face, the flux is formulated as 

f = 5 < n > 

where the dissipative flux, d, is constructed from a 
blending of a first-difference and pseudo third-differ- 
ence. When integrated over the faces of the control vol- 
ume, this yields a blended Laplacian and Biharmonic 
operator which is used to dissipate spurious oscillations. 
This scalar dissipative flux is formulated as 

d = d l -d . 3 

d { = 8e< 2 > {W R -W L ) (12) 

d 3 = ae< 4 ) [D 2 (W r )-D 2 (W l )] 

The pseudo-second differences formed for each cell, 


5 

American Institute of Aeronautics and Astronautics 



D 2 , represent undivided approximations to Laplacians 
using only differences about faces. That is, for an arbi- 
trary function, f 

D 2 (f) = X (f R -f L ) (!3) 

faces 

where the R and L refer to the logically oriented cells 
sharing the face. 

The coefficients £ ^ and £ ^ are used to scale the 
first and pseudo-third differences with the cell size, and 
to turn off the fourth-order dissipation across shocks. 
The coefficient of the first-difference in (12) is com- 
puted as a normalized second-difference of pressure, 
acting as a shock sensor 


D^jP) 

I Vr+'i> 

faces 


This gives the coefficients 


/ 2 ^ = K 2 max(v^,v L ) 
£< 4) = max[0,(K 4 -£ (2) )] 


(14) 


(15) 


Standard values of the dissipation coefficients are used 
as = 1/2 and k 4 = 1/32 .The factor c scales 
the dissipation according to the maximum eigenvalue of 
the flux jacobian as 


(cj+a*) 

2 

°l/r= (i“V v v w/ y +a >i L/ * 


(16) 


On a uniform mesh, this dissipation results in a blended 
Laplacian/biharmonic operator. 

IILd Treatment ofViscous Fluxes 

The viscous fluxes are computed following a directed 
gradients procedure suggested by Mitchell [23]. This 
procedure is linearly K-exact, and produces the same 
gradient computed using a divergence-based reconstruc- 
tion that preserves linear data, as in [24] and [4], com- 
monly called a diamond-path reconstruction. The idea is 
based upon taking the inner product of the gradient with 
three vectors joining locations where the data has been 
obtained in an at least linearity-preserving manner. That 
is, for the three vectors 8 j, 8 2 , 83 


A = (Vk) • 8j 

A« 2 = (Vk) • 8 2 (17) 

A « 3 = (Vw) ■ 8 3 

This results in a 3x3 linear system for the unknown gra- 
dient, as 


Axj Ay l Az x 


u 

X 


A Mj 

Ax 2 Ay 2 A z 2 


u y 

= 

A “ 2 

Ax 3 Ay 3 Az 3 


“z 


A m 3 


The choice of the base vectors, 8 f is taken so that 8 ^ 
joins the two centroids of the cells that share the face, 
and the two others lie in the plane of the face. When a 
linearly-exact procedure is used to produce the data at 
the vertices, this procedure preserves linear gradients. 
The viscous flux construction, evaluation and scattering 
to the cell residuals is done on a face basis, which for 
vectorization, depends upon loop coloring, as is 
explained in the next section. 

IV Loop Coloring for Vectorization 

Vectorization of the solver is pursued for the use of 
CRAY class vector machines. If the solver is used upon 
non-vector machines, such as workstations, or is made 
parallel, the basic structure and loop indexing used for 
parallelization will result in a negligible penalty. The 
preprocessing time for the loop colorings is marginal. 

There are primarily two different types of scatter opera- 
tions that must use loop coloring to obtain vectorization. 
Loops over faces scattering to cells are performed for 
the inviscid reconstruction, time-step computation and 
for the viscous and inviscid flux integrals. Loops over 
cells, scattering to vertices are performed to interpolate 
the data at the vertices from the cells. All loops which 
are to be vectorized are done so by compiler directives 
to ignore vector dependencies. The vectorized code 
ports to non- vector machines with no changes. 

The face loop coloring is performed in a heuristic fash- 
ion as indicated by the following pseudo-code. 

AT_ color s**0 

/* loop over all faces */ 
for (all faces) 

added_color=fal se 

/* loop over all available colors */ 
for (all colors && ! added_color) 

L_clrd=is_left_colored_color ( color) 
clrd=± s_rigb t_co 1 ore d_col or ( color) 

/* if left and right cell ore not color 
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color them color */ 
if( !L_clrd ££ lR_clrd) 
face_color=color 
edded_color=true 
endif 
endfor 

/* if not colored, color it color */ 
if ( l added., col or) 

face_color=N_col ors 
N_ colors** 1 
endif 
endfor 

This results in roughly even-length loops of faces, 
except for the last one or two colors, which typically 
have a smaller numbers of faces. The balance of the col- 
ors is dependent upon the order that the above algorithm 
visits the faces, but using the input ordering appears to 
be sufficient. 

Vectorization of the cell-to-vertex scattering operation is 
complicated by the fact that the cells in the mesh may 
have a variable number of vertices. This is overcome by 
grouping cells into groups of like vertex number, termed 
here as verticity. Within cells of the same verticity, cells 
are colored to avoid scattering to the same vertex in the 
same loop. A binary-encoded coloring history of each 
vertex is used to group the cells into colors. The follow- 
ing pseudo-code illustrates the approach. 


/* zero history for ell vertices */ 
for (ell vertices) 
hlst**0 
endfor 

/+ loop over like-vertici ty cells */ 

nColors*0 

for (ell cells ) 

/* construct e single history for this 
cells' vertices using bitwise or V 
v_C*Q 

for (ell vertices in cell) 
v_C*vjC I hist 
endfor 

f_a_c=first_evei l^color ( v_C) 
num*l 

for (ell colors < f_e_c) 
num***2 
endfor 

for (ell vertices in cell) 
hist*hist I num 
endfor 
endfor 

The routine first_eveii_coior performs bitwise shifts 
upon the cell-encoded vertex history and finds the first 
available color for the vertices in the cell. 


int first_evel_color(vC_hist) 


POB m X « 0 
colors 0 

while ( (vC_hist £ pos) 0) 
pos=pos « 1 
color** 
endwhile 
return color 

This colors the loops within a verticity-group of roughly 
even length loops, with the last one or two loops of 
smaller length. Representative performance values of 
the vectorized code are given in the following sections. 


V Validation and Demonstration 

The accuracy of the solution procedure is first assessed 
by performing a grid refinement study upon a sequence 
of related meshes of hexahedra, prisms, tetrahedra and 
those produced by a Cartesian grid generator. A devel- 
oping laminar boundary layer and a developing laminar 
duct flow are shown. Also shown is the laminar flow 
computed about a flat plate/cylinder intersection, upon a 
prismatic/tetrahedral grid generated by Connell [8]. 

V.a Grid Refinement Study: Supersonic. 
Vortex 

An analytical solution to the Euler equations is used to 
assess the accuracy of the convective flux discretization 
schemes of the solver. This flow field has been used by 
Aftosmis et al. [25] to investigate the effects of gradient 
limiting and by Luo, et al. [26] to assess an improved 
reconstruction scheme for triangular meshes. The solu- 
tion is a relatively simple function, and also lends itself 
well to a grid refinement study, since although the flow 
is simple, there are considerable gradients in pressure 
and density across the domain. By assuming an isen- 
tropic flow described by a line vortex situated at the ori- 
gin, parallel to the z-axis, the solution is 




= 1 + 


p _ 

\ ' 

/ry/(Y-i) 

Pi 

U.J 

p 

fT\ Y/(Y-1) 

p i = 

U ( .J 


u = U cos0 
v = UsinQ 
w = 0 


( 19 ) 


where the i-subscript refers to conditions along a refer- 
ence (inner) radius r. , % = r/r. the polar angle in the 
z=constant plane is 0 = atan {y/ x) , and the flow 
speed non-dimensionalized by the inner radius speed is 
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U = l/£ . The domain is taken to be for 
[1, 1.384], z€ [0,1] and 0€ [0 ,tc/ 2] , and 
the Mach number along the inner radius is taken to be 
= 2.25. This is a two-dimensional problem, so dur- 
ing the refinement sequence, the discretization of the 
grid in the z-direction is kept the same for all of the 
grids. Figure 1 shows Mach number contours for this 
flow. 



Figure 1. Mach contours, compressible 
vortex flow. 

A sequence of related grids are generated based upon a 
sequence of structured hexahedral grids with 5x5x5, 
10x10x5, 20x20x5 and 40x40x5 cells. Tetrahedral grids 
are created from the hexahedral meshes by creating six 
tetrahedra from each hexahedra. Three families of pris- 
matic grids are generated from the base hexahedral 
meshes by creating prisms from the hexahedra whose 
orientation of the normal vector of the triangular faces 
in the prisms lie along the three different computation 
coordinate axes. A sequence of Cartesian grids are also 
generated, which are also used to assess the accuracy, 
but are not directly related to the five grid sequences 
above. Figure 2 shows the relationships between the 
hex-related meshes, while Figure 3 shows the intersec- 
tion of a z=constant plane through a representative Car- 
tesian mesh. Note that for plotting purposes, cut- 
Cartesian cells that cannot be represented as either hexa- 
hedra, pyramids, prisms or tetrahedra are split into a col- 
lection of tetrahedra and pyramids. For these cut cells, a 
fictitious point located at the cell centroid is introduced, 
which is used to create, corresponding to each triangular 
or quadrilateral face of the cell, a tetrahedron or a pyra- 
mid, respectively. This procedure is only needed for 
plotting purposes, and is not used in the flow solver. 




Prisms, p_id=2 Prisms, p Jd=3 



(Base) Hexahedron 


Figure 2. The five, hexahedra-derived grids 
types used in the mesh convergence study. 



Figure 3. z=constant cut through a 
Cartesian mesh 


A plot of the Lj norm of the density error 
e = |P ” Pexacrl a g a ^ nst a representative two- 

dimensional length scale found on each mesh is shown 
in Figure 4. This length scale is introduced merely for 
the estimation of the truncation error order, and is not 
representative of a computational cost or efficiency 
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norm. The intent of this study is only to gauge the accu- 
racy and correctness of the inviscid flux constructions 
and is not an attempt to assess whether one mesh topol- 
ogy is superior to another. That assessment would 
require a careful comparison of cost as well. This repre- 
sentative length scale is calculated as 


l 


2d ~ 



y nCells_ 

ave ~ nCells 


( 20 ) 


All calculations were performed using the upwind for- 
mulation Ill.b, with the FDS scheme without gradient 
limiting. 



The orders of the discrete truncation error, found by 
computing the slope on a logarithmic plot of the final 
two meshes in each sequence, for the six different mesh 
sequences, are shown in Table I, in the Appendix. All 
schemes were globally second-order accurate, while the 
tetrahedral and plane_id=l meshes had a first-order 
max-norm. 

When viewing Figure 4, it is important to keep in mind 
that the hexahedral, tetrahedral and prismatic meshes 


are all interrelated, since they are derived from the same 
mesh, while the Cartesian is not. All meshes are made 
with a constant spacing in the z-direction, resulting in a 
stack of four cells in the z-direction. The Cartesian mesh 
and the hexahedral mesh are shown to have nearly the 
same truncation error for the same length scale, which is 
attributable to the similarity in the reconstruction/flux 
construction schemes that they use. In [27] an analytic 
solution of the Euler equations, Ringleb’s flow, was 
used to compare the error computed by the Cartesian 
approach and a structured mesh approach. There, a 
structured mesh, which used an upwind coordinate-by- 
coordinate reconstruction, was shown to be slightly 
more accurate than the Cartesian approach, which used 
a multi-dimensional reconstruction. The results shown 
here compare the truncation error between a sequence of 
hexahedral and Cartesian meshes, where the multi- 
dimensional reconstruction shown in Section IIIb.1 is 
used. 



A series of calculations were also made using the cen- 
tral-difference flux scheme with added dissipation on 
the exact same meshes. Figure 5 shows the computed 
L norms against the two-dimensional length scale and 
Table II shows the asymptotic orders of accuracy. 

As is indicated by the order of the L ^ norms, the partic- 
ular implementation of the dissipative flux constructions 
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is not uniformly second-order accurate. It is surmised 
that the construction of the dissipative fluxes near the 
boundaries has reduced the accuracy of the scheme. The 
results shown here construct the third-difference opera- 
tor at the boundaries using ghost-cell information. A 
comparison to a presumably less accurate formulation 
which does not use ghost cell data indicated no appre- 
ciable improvement. 

Since the discrete errors for both the upwind and central 
difference schemes are available for the same meshes, 
one can also compare these two schemes. Figure 6 
shows the L ^ norms for the hexahedral, tetrahedral and 
Cartesian meshes. 



Figure 6. Comparison of L x -norms for 
the upwind and central schemes. 


As is seen in the plot, there is an appreciable difference 
between the truncation error of the central and upwind 
schemes on the Cartesian mesh, where the central 
scheme exhibits nearly ten times the error on the finer 
mesh than the upwind scheme. This difference in error 
is not at as severe on the tetrahedral and hexahedral 
meshes, where the central scheme error is approxi- 
mately four and six times more than the upwind formu- 
lation, respectively. 

The grid convergence study was performed on IBM 
RS6000, model 590 workstations, using the xlc com- 
piler with standard optimizations. The processing rates 


for both the upwind and central schemes are shown in 
Table III and Table IV for the final meshes in each 
sequence. Also shown in are the rates for the same 
meshes and algorithms, but on the CRAY C-90, eagle, 
using the Cray standard C compiler. Version 4.O.2.7. 

It is important to consider the particular way the differ- 
ent schemes are compared to one another in Figure 4, 
Figure 5 and Figure 6. The length scale has been con- 
structed to deduce the order of the schemes. A compari- 
son that would critically gauge the different approaches 
upon the different mesh topologies must use some mea- 
sure of computational efficiency, perhaps measured by a 
memory-time integral. This efficiency will be greatly 
effected by the tailoring of the solver to a particular con- 
vective flux discretization (i.e. central or upwind) and 
for a particular mesh topology. The intent of this conver- 
gence study is only to verify the accuracy of the particu- 
lar flux constructions and to validate the mixed volume 
grid approach, and not to promote one mesh topology 
over another. 

V.b Developing Laminar Flow over a Flat. 
Plate 

The flow developing over a flat plate which is aligned 
with the free stream is computed next A mesh of hexa- 
hedra is generated with 31, 21 and 5 points in the 
streamwise, plate normal, and span- wise directions, 
respectively. The mesh is stretched in both the stream- 
wise and normal directions. Flow conditions correspond 
to a freestream Mach number of = 0.25 and a 
Reynolds number based on plate length of 
Re = 50, 000 . The code is ran using the upwinding 
formulation and the FDS scheme. As in the previous 
case, a hexahedral mesh has been constructed, and from 
this mesh, other derived mesh types are found. A tetra- 
hedral mesh and a prismatic mesh with the p_id= 2 are 
constructed from the base mesh, and the prismatic mesh 
is shown in Figure 7. 



Figure 7. Prism (pjd=2) mesh derived 
from hexahedral mesh. 
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Figure 8 and Figure 9 show the computed u- and v- 
velocity profiles plotted against the similarity coordinate 
r\ , for the three meshes at a location on the mid-plane of 
the plate, with a Reynolds number based on distance 
from the leading edge of Re x = 40, 000 . Each mesh 
produces good results. 



U/U 

oo 

Figure 8. u- velocity profiles at 

Re x = 40, 000 . 


V.c Developing Laminar Flow in a 
Rectangular Duct 

The developing laminar flow in a five-to-one aspect 
ratio duct is computed next, and compared to the experi- 
mental data of Sparrow et al. [28]. The computations are 
performed for a Reynolds number based on inlet mass 
flow rate and hydraulic diameter of 500 for a duct length 
of 31.25 inches. The 5:1 aspect ratio duct has major and 
min or widths of 3.125 and 0.625 inches. The exit pres- 
sure is specified according to the experimental data from 
the experimentally measured pressure drop correlation 


K 




1.89 + 76.3 



( 21 ) 


where D is the hydraulic diameter. The inlet total condi- 
tions are set so that the Mach number based on area- 
averaged axial velocity is M m = 0.1 . A 15x15x21 
hexahedral mesh is constructed with stretching in all 
directions. A tetrahedral mesh is also constructed, as 
before, from the base hexahedral mesh. These calcula- 
tions use the upwinding formulation with the FDS flux 
function. Figure 10 and Figure 11 compare the com- 
puted centerline pressure drop and streamwise velocities 
from the hexahedral and tetrahedral meshes to the 
experimental data of [28]. Agreement with experiment 
is considered good. 



Figure 9. v-velocity profiles at 
Re x = 40, 000 



Figure 10. Computed and experimental 
axial-pressure loss. 
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V.d Wall/Cvlinder Flow: Prismatic/ 
Tetrahedral Mesh 

The flow about a wall/cylinder combination is computed 
using the mixed volume grid approach. The mesh was 
generated by the prismatic/tetrahedral mesh generator 
developed by Connell et al. [8], for which the solver 
provides an input data type. Flow conditions correspond 
to a Reynolds number based upon cylinder diameter of 
Re = 50 and a freestream Mach number of 
M = 0.25 . Figure 12 shows a perspective view of the 
grid, while Figure 13 shows a slice through the mesh at 
the mid-plane, showing the prismatic mesh about the 
cylinder. The mesh has 9810 prisms and 19243 tetrahe- 
dra. The calculation converged approximately four 
orders of magnitude in 10,000 iterations, using 5011 
seconds of cpu time and approximately 4.8 Mwords of 
storage on the Cray C-90, eagle. Figure 14 shows a per- 
spective view of speed contours of the computed solu- 
tion at various planes slicing the volume mesh, showing 
the boundary layer growth along the wall, the upstream 
influence of the cylinder and the momentum deficit in 
the wake. The predicted solution appears to be symmet- 
ric. 



Figure 12. Perspective view of 
boundaries of prismatic/tetrahedral 
mesh for wall/cylinder. 



Figure 13. Upper boundary of prismatic/ 
tetrahedral mesh 



Figure 14. Speed contours through 
volume 
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VI Concluding Remarks and Future Efforts 

An approach which solves the compressible Euler and 
Navier-Stokes equations upon control volumes of nearly 
arbitrary polyhedra has been presented. The conserva- 
tion laws are solved using a cell-centered, finite-volume 
formulation where the control volumes are created from 
an arbitrary number of triangular and quadrilateral 
faces. This flexibility allows the unified treatment of 
structured, unstructured, prismatic, and Cartesian-cell 
based grids with a single flow solver. The basics behind 
the organization and interrogation of the data structures 
has been explained. 

The mixed volume grid approach uses both cell-aver- 
aged and vertex data. The convective terms of the gov- 
erning equations have been treated with both an upwind 
and central differencing approach. The reconstruction 
method used for the upwind flux discretization uses an 
application of the divergence theorem upon the actual 
control volume itself. The viscous fluxes have been con- 
structed using a directed gradients procedure. 

An analytical solution to the Euler equations, a super- 
sonic vortex, has been used to assess the accuracy of the 
approach, and has been used to compare results amongst 
the different grid types and between the upwind and 
central difference approximations. A refinement 
sequence of four hexahedral meshes was used to gener- 
ate 3 families of prismatic meshes and a sequence of tet- 
rahedral meshes by subdividing each hexahedra into the 
respective volume type. In addition, a Cartesian gener- 
ated mesh was used and compared to the others. Here, 
the bulk of the Cartesian generated volume is comprised 
of axes-aligned hexahedra with arbitrary polyhedra 
along the boundaries created from the geometric inter- 
section of the boundaries and the Cartesian cells. The 
grid refinement study showed that the hexahedral based 
mesh had the lowest error, while the tetrahedral mesh 
had the highest In addition, the global error norms indi- 
cated that for the same mesh the upwind formulation 
was more accurate than the central difference approach. 
The Cartesian grids were shown to give an error compa- 
rable to the hexahedral meshes. The upwind formulation 
was shown to be globally second-order accurate upon 
all the meshes, and locally second-order accurate on 
some of the element types. The central scheme was 
shown to be globally, marginally second-order accurate, 
and first-order accurate locally. This reduction in the 
order for the central scheme was surmised to be caused 
by the construction of the boundary dissipative fluxes. 

Laminar calculations have been made for the develop- 
ing flow over a flat plate for hexahedral, prismatic, and 
tetrahedral meshes and for the developing flow in a 5:1 
aspect ratio duct. Favorable comparisons to theory and 


experiment were shown for these cases. To demonstrate 
a prismatic mesh type calculation, a prismatic/tetrahe- 
dral mesh, supplied by Connell et al. [8], was used to 
compute the flow about a cylinder/wall configuration. 
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A ppendix A: Mesh Convergence Study: 
Asymptotic Slopes 

TABLE I Asymptotic Slopes of Error Norms: 


Upwind Formulation 


Grid Type 

L 1 

Norm 

L 

oo 

Norm 

Hexahedral 

2.08 

1.94 

Tetrahedral 

1.77 

0.93 

Prismatic : plane_id= 1 

1.68 

0.88 

Prismatic: plane_id=2 

2.05 

1.63 

Prismatic: plane_id=3 

2.02 

1.78 

Cartesian: Az = 1/4 

2.03 

1.73 


TABLE II Asymptotic Slopes of Error Norms: 


Central Differencing Formulation 


Grid Type 

L i 

Norm 

L 

oo 

Norm 

Hexahedral 

1.82 

1.18 

Tetrahedral 

1.17 

0.78 

Prismatic: plane_id=l 

1.85 

0.78 

Prismatic: plane_id=2 

1.33 

1.02 

Prismatic: plane_id=3 

1.82 

1.02 

Cartesian: Az = 1/4 

1.96 

1.06 
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A ppendix B; P rocessing Rates on FinaL 
jUeshes used in Conve rgence StudY 


TABLE m Cell processing rates (seconds/ 
iteration/cell) on final meshes, upwind scheme 


Grid Type 

RS6000 

Model 

590 

CRAY 

C-90, 

eagle 

Hexahedral 

3.41e-4 

2.51e-5 

Tetrahedral 

2.27e-4 

1.50e-5 

Prismatic: plane_id=l 

2.88e-4 

2.04e-5 

Prismatic: plane_id=2 

2.89e-4 

2.02e-5 

Prismatic: plane_id=3 

2.79e-4 

2.05e-5 

Cartesian: A z =1/4- 

3.47e-4 

2.20e-5 


TABLE IV Cell processing rates (seconds/ 
iteration/cell) on final meshes, central scheme 


Grid Type 

RS6000 

Model 

590 

CRAY 

C-90, 

eagle 

Hexahedral 

1.84e-4 

1.94e-5 

Tetrahedral 

1.26e-4 

1.14e-5 

Prismatic: plane_id=l 

1.55e-4 

1.57e-5 

Prismatic: plane_id=2 

l.53e-4 

1.54e-5 

Prismatic: plane_id=3 

1.50e-4 

1.55e-5 

Cartesian: Az = 1/4 

1.82e-4 

1.60e-5 
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